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STORAGE  FLOOD  ROUTING  WITHOUT  COEFFICIENTS 


D.  L.  Brakensielc^ 


INTRODUCTION 

Flood  routing  methods  that  are  based  on  the  conseirvation  of  mass  (continuity)  and  a  single- 
valued  relationship  between  depth  of  flow  and  flow  rate  are  usually  termed  "storage  flood 
routing"  methods.  This  type  of  formulation  is  discussed  more  rigorously  under  the  theory  of 
kinematic  waves.  A  fundamental  reference  to  this  theory  Is  made  by  Lighthill  and  Whitham  (7). 
Henderson  (4)  considered  the  kinematic  treatment  of  flood  waves  in  prismatic  channels.  More 
recently,  Henderson  and  Wooding  (,5)  applied  the  kinematic  wave  theory  to  the  overland  flow 
problem. 

Many  hydrologists  using  storage  flood  routing  utilize  coefficients,  such  as  X  and  K  of  the 
Muskingum  method  (2)  or  C  of  the  convex  method.  This  report  presents  a  numerical  method 
for  storage  flood  routing  without  coefficients.  A  number  of  flow  problems  are  solved  with  a 
computer  program.  Some  limitations  of  the  kinematic  formulation  are  Illustrated  by  example. 

FORMULATION 

The  formulation  of  storage  of  kinematic  flood  routing  is  composed  of  the  continuity  equation 


6Q      ,.    ^A 
3x  ^t 


=    q  (1) 


and  a  rating  function 

Q    =    Q(A)  (2) 


where 


Q    =    flow  rate,  c.f.s., 

2 

A    =    flow  area,  ft.  , 

q     =    lateral  Inflow  (+)  or  lateral  outflow  (-),  c.f.s./ft.,  and 
X,  t     =    distance  (feet)  and  time  (seconds)  coordinates. 


*  Research  hydraulic  engineer.  Soil  and  Water  Conservation  Research  Division,  Agricultural  Research  Service, 
U.S.  Department  of  Agriculture. 

'Underscored  numbers  in  parentheses  refer  to  Literature  Cited,  page  20. 

•Mockus,  Victor,  Chapter  II,  Computer  program  for  project  formulation  hydrology.  Prepared  for  Soil  Conser- 
vation Service,  U,S,  Department  of  Agriculture,  by  C-E-I-R,  Inc.  January  1964, 


The  salient  features  of  the  continuity  equations  are  obtained  by  substituting  Q  =  AV  into  equa- 
tion   1,    differentiating    the    product,    and   multiplying   both   sides  by  At  Ax,  to  obtain  volumes. 


(a) 


(b) 


AxAt       -^:: +    -= + 


ax 


^x 


(c) 
at 


=     q  (AxAt). 


where 

V  =  velocity,  ft./sec. 
A  t  =  time  increment,  sec. 
Ax  =  length  increments,  ft. 

The  terms  on  the  left-hand  side  are  usually  referred  to  as  follows: 

(a)  -  prism  storage, 

(b)  -  wedge  storage,  and 

(c)  -  time  rate  of  change  of  storage. 


The  usual  treatment  of  equations  1  and  2  is  to  recast  them  as  follows: 
Using  the  grid  pattern  in  figure  1,  equation  1  can  be  written  (wi±  q  =  O), 
± 


Lt 


At 


At 


Ax  Ax  Ax 


Figure  l.~Grid  for  finite  difference  approximations. 


->-  X 


dQ  Q4-  Qt 

dx  "  Ax 

<3A  _  A4  -  A3  f  Ag  -  A| 

at  ""  2At 


q    = 


q4  +  q2 


Q4  -    Q2    +    (A4  +   A2)      -      (A3    +     Aj )      ^     q 
Ax  2At 


or 


Now  if  we  define 

S2    =(^^    ^    ^^^  )ax     =    storage  at  t   =    t    +  A  t, 

Sj^    =-L_ 3_^ lL  Ax    =  storage  at  t  =  t, 

Q_^  =  O  =  outflow,  and 
Q2  =  I    =  inflow. 

The  equation  3a  can  be  written 

S2     -     Si 


(3a) 


I    -    O    = 


At 


or 

(I    -     O)  At    =   AS.  (3b) 

Equation  3b  is  similar  to  the  usual  continuity  expression  (inflow  -  outflow  =  change  of  storage). 
Storage  is  evaluated  as 

S    =     S(0)  (4) 

or 

S    =     S(0,I)  (5) 

where 

S    =     storage,  c.f.s.-hrs., 
O    =     outflow,  c.f.s.,  and 
I    =     inflow,  c.f.s. 

Equation  4  is  used  by  methods  that  neglect  wedge  storage,  such  as  the  storage-indication  method. 
Equation  5  is  used  by  the  Muskingum  method  where  two  coefficients  X  and  K  are  introduced  to 
account  for  both  prism  and  wedge  storage.  The  development  in  this  report  is  based  on  only 
equations  1  and  2,  i.e.,  no  functional  form  is  assumed  for  relating  reach  storage  to  inflow  and 
outflow.  The  assumption  inherent  in  kinematic  wave  theory  is  not  altered;  i.e.,  the  equation  of 
motion  can  be  approximated  by  a  single-valued  function,  equation  2. 


SOLUTION  FORMULATION 

A  simultaneous  solution  of  equations  1  and  2  is  accomplished  by  obtaining  a  numerical 
solution  on  a  small  computer  (IBM  1620,  40  K  memory  and  indirect  addressing).  Equation  1  is 
approximated  by  the  finite  difference  quotients  shown  in  figure  1.  Substituting  these  quotients 
into  equation  1,  one  obtains: 

Q4    -     Q2      ,     A4     -     A3     +     A2     -     Ai     _     - 
— ^^ "^  2Kt  -    ^ 


or 


-ALq,     .41-=    (^^^),    A^Q,.    (^^A.,) 


or  ,  , 

XQ,    +     4^    =     a+  P  (6) 


where 

X  =    At/Ax,  .      ;•- 

a   =     (Ai     +  A3)   /2,  and 
fi   =     XQ2    +    /-A2    +    Atq\ 

The  subscripts  in  these  equations  are  defined  as  follows: 

A-|^  =  flow  area  at  location  x  and  time  t, 

A3  =  flow  area  at  location  x  +Ax  and  time  t, 

Ag,  Q2  =  flow  area  or  rate  at  location  x  and  time  t  +  At,  and 

A4,  Q4  =  flow  area  or  rate  at  location  x  +  Ax  and  time  t  +  At. 

Equation  2  is  retained  in  the  computer  memory  as  an  array  of  numerical  values,  i.e.,  paired 
Q  and  A  values. 

The  rating  functions  are  determined  by  one  of  several  procedures.  If  normal  flow  (turbulent 
flow)  is  assumed,  the  Manning  or  Chezy  equation  can  be  applied  at  each  channel  section.  An 
alternative  for  natural  channels  is  to  calculate  a  series  of  water  surface  profiles;  e.g.,  using 
U.S.  Department  of  Agriculture  Hydrograph  Laboratory  Computer  Program  No.  2  (8).  The  com- 
puted values  for  rates  of  flow  and  depths  of  flow  at  each  section  along  the  channel  system  define 
a  section  rating  function. 

Since  the  rating  function  is  not  an  explicit  function  but  a  set  of  numerical  values,  a  simul- 
taneous solution  of  equations  6  and  2  requires  an  interation  procedure.  In  this  study,  the  pro- 
cedure used  is  known  as  the  method  of  false  position  6.  A  special  constraint  is  required  to  account 
for  incipient  flow  at  downstream  sections  during  rising  stages  and  for  cessation  of  flow  at  up- 
stream sections  during  falling  stages,  i.e.,  for  value  of  (a  +  /S  )  <  O,  Q4  and  A4  are  set  equal 
to  base  flow. 


SOLUTION  PROGRAM 

Presented  in  appendix  1  are  the  source  programs  for  the  main  line  and  three  subroutines. 
These  programs  are  written  in  Fortran  II. 

Program  input  and  output  formats  are  presented  in  appendix  2.  Arrays  of  data  are  read  in 
by  way  of  the  "Read"  subroutine.  Flow  area  values  are  calculated  for  flow  rate  values  by  an 
interpolation  function.  Arguments  are  obtained  from  the  "Table  Look-up"  subroutine.  One  limita- 
tion on  inflow  input  is  that  the  final  rate  of  inflow  must  equal  the  initial  value  (return  to  a  dry 
channel  or  to  the  initial  base  flow).  Note  that  lateral  flow  hydrographs  must  have  the  same 
number  of  entries  as  the  inflow  hydrograph.  Dummy  zeros  can  be  entered. 

In  storage  flood  routing  methods,  calculations  proceed,  as  usual,  in  a  downstream  direction. 
An  inflow  hydrograph  is  given  at  the  upstream  section  as  a  boundary  condition.  For  overland  flow 
applications,  this  is  the  "null"  hydrograph.  The  initial  conditions  are  known  rates  of  flow  and 
flow  areas.  These  can  be  zero  if  a  dry  channel  is  assumed.  If  base  flow  is  assumed,  flow  area 
should  be  calculated  from  the  section  rating  functions.  Lateral  flow  (q)  can  be  put  in  as  a  plus 
quantity  for  local  inflows  or  as  a  minus  quantity  for  transmission  losses.  Inflow  from  a  tributary 
can  be  treated  as  lateral  flow  introduced  over  a  short  channel  length. 

Initially  the  program  section,  "Routing  During  Inflow,"  determines  the  time  of  incipent  flow 
at  the  next  downstream  section;  thereafter,  inflow  is  routed  by  application  of  the  "Solve"  sub- 
routine. This  subroutine  performs  the  iteration  solution  of  equation  6  with  discharge  values  cal- 
culated from  the  rating  function  table.  Iteration  cycling  is  terminated  by  a  tolerance  (TOLR)  level. 

After  the  last  value  of  inflow  has  been  routed,  computations  are  transferred  to  the  main 
line  program  section  "Routing  After  Inflow"  section  of  the  program.  Routing  computations  con- 
tinue until  downstream  outflow  decreases  to  a  value  within  0.1  of  the  initial  value;  the  last  out- 
flow entry  is  then  set  equal  to  the  initial  inflow  value.  All  input  data  arrays  are  filled  out  to 
correspond  in  length  to  the  outflow  arrays.  During  solution  punch  out,  outflow  section  data  are 
interchanged  with  inflow  section  data  so  as  to  become  inflow  data  for  the  next  routing  reach. 
After  data  arrays  are  read  in  for  the  next  outflow  section,  computations  proceed  for  the  next 
routing  reach.  Problem  solution  is  complete  after  the  last  section  data  have  been  read  in  and 
a  Reader  "no  feed"  signal  occurs  on  the  IBM  1620  console. 

EXAMPLES 

Example  1  is  an  inflow  hydrograph  introduced  on  a  constant  base  flow  of  150  c.f.s.  Manning's 
equation  was  utilized  for  the  rating  function  of  a  theoretical  rectangular  channel  30  feet  wide, 
with  a  bed  slope  of  0.001  and  an  n-value  of  0.03.  Since  this  was  a  prismatic  channel,  the  rating 
functions  are  the  same  at  all  sections  of  the  channel.  The  time  scale  increment  (At)  was  3,600 
seconds,  and  the  distance  scale  increment  (Ax)  was  7,200  feet.  Table  1  presents  program  input 
and  table  2  presents  the  program  output.  Figure  2  is  a  succession  of  routed  hydrographs.  Figure 
3  is  a  typical  storage  outflow  function  generated  by  the  routing  method.  The  storage  values  were 
calculated  by  the  average  end- area  method.  The  storage  loop  indicates  the  importance  of  wedge 
storage.  If  the  wedge  storage  term  in  the  continuity  equation  were  neglected,  the  relationship 
would  be  a  single-valued  curve. 

For  example  2,  an  inflow  hydrograph  was  introduced  on  an  initially  dry  channel  with  a 
constant  (in  time  and  space)  lateral  outflow  (q)  of  0,05  c.f.s./ft.  The  same  channel  geometry  was 
used  as  in  example  1,  and  Manning's  equation  was  used  to  calculate  the  rating  functions.  The 
program  input  Is  presented  in  table  3.  A  succession  of  routed  hydrographs  is  presented  in 
figure  4.  The  rather  uniform  decline  in  the  hydrograph  results  from  the  constant  transmission  loss. 

Figure  5  presents  two  additional  surface  flow  problems.  Presented  in  figure  5  at  top,  is  a 
set  of  hydrographs  resulting  from  a  constant  inflow  of  150  c.f.s.  This  problem  would  be  typical 
of  the  border  of  furrow  irrigation  problem.  Also,  releases  from  a  detention  reservoir  might 
lead  to  a  similar  problem.  In  this  example,  transmission  losses  were  set  to  zero.  The  program 
input  for   this   example   is   given   in   table  4.  Figure  5,  at  bottom,  illustrates  the  overland  flow 


Table  1. — Input  for  example  1 
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Figure  2. — Succession  of  routed  hydrographs. 


Table  2. — Output  for  example  1 


IN  SECTION  NO  » 

1 

OUT  SECTION  NO 

=   2 

IN  AREA 

IN  DISCH 

OUT  AREA 

OUT  DISCH 

62.915 

150.000 

62.915 

150.000 

118.357 

400.000 

99.849 

307.579 

162.843 

650.000 

147.341 

557.596 

231.783 

1088.000 

211.841 

954.559 

328.012 

1775.000 

302,998 

1587.613 

456.162 

2775.000 

424.965 

2524.881 

576.877 

3775.000 

547.966 

3531.256 

635.632 

4275.000 

620.306 

4143.905 

650.242 

4400.000 

645.582 

4360.114 

635.632 

4275.000 

638.202 

4296. 9ft8 

584.340 

3038.000 

595.476 

3932.010 

524.637 

3338.000 

538.f.97 

^454.469 

456.162 

2775.000 

472.599 

2909.573 

401.534 

2338.000 

415.538 

2449.689 

353.016 

1963.000 

365.748 

2061.307 

311.326 

1650.000 

322.645 

1734. 7B9 

268.029 

1338.000 

280.105 

1423.800 

231.783 

1088.000 

242.763 

1161.472 

203.668 

900.000 

212.420 

958.435 

173.232 

713.000 

183.072 

772.670 

141.375 

525.000 

152.521 

587.407 

118.357 

400.000 

127.306 

448.132 

92.637 

275.000 

103.261 

324.558 

62.915 

150.000 

76.982 

205.916 

62.915 

150.000 

65.720 

161.147 

62.915 

150.000 

63.475 

152.223 

62.915 

150.000 

63.027 

150.443 

62.915 

150.000 

62.915 

150.000 

IN  SECTION  NO  « 

2 

OUT  SECTION  NO 

=   3 

IN  AREA 

IN  DISCH 

OUT  AREA 

OUT  DISCH 

62.915 

150.000 

62.915 

150.000 

99.849 

307.579 

86.461 

247.099 

147.341 

557.596 

130.575 

465.990 

211.841 

954.559 

192.227 

828.188 

302.998 

1587.613 

278.227 

1410.453 

424.965 

2524.881 

394.877 

2286.205 

547.966 

3531.256 

518.228 

3284.904 

620.306 

4143.905 

602. 050 

3987.742 

645.582 

4360.114 

638.380 

4298.508 

638.202 

4296.988 

638.993 

4303.754 

595.476 

3932.010 

604.556 

4009.174 

538.697 

3454.469 

551.875 

3563.929 

472.599 

2909.573 

488.422- 

3039.118 

415.538 

2449.689 

430.014 

2565.149 

365.748 

2061.307 

378.825 

2162.272 

322.645 

1734.780 

334.335 

1822. "^64 

Table  2. — Output  for  example  1 — Continued 


280.105 

1423.800 

292.033 

1508.545 

242.763 

1161.472 

253.818 

1237.030 

212.420 

958.435 

221.745 

1020.835 

183.072 

772.670 

192.665 

830.846 

152.521 

587.407 

162.528 

648.094 

127.306 

448.132 

136.630 

499.075 

103.261 

324.558 

112.799 

372.226 

76.982 

205.916 

88.473 

256.190 

65.720 

161.147 

72.484 

188.034 

63.475 

152.223 

65.718 

161.139 

63.027 

150.443 

63.653 

152.931 

62.915 

150.000 

63.085 

150.673 

62.915 

150.000 

62.949 

150.134 

62.915 

150.000 

62.915 

150.000 

IN  SECTION  NO  ■ 

3 

OUT  SECTION  NO 

=   4 

IN  AREA 

IN  DISCH 

OUT  AREA 

OUT  DISCH 

62.915 

150.000 

62.915 

150.000 

86.461 

247.099 

77, 

.694 

208, 

.744 

130.575 

465.990 

115, 

►  222 

384, 

►  332 

192.227 

828.188 

172, 

.598 

709, 

.159 

278.227 

1410.453 

254, 

.583 

1242, 

.469 

394.877 

2286.205 

365, 

.414 

2058, 

.725 

518.228 

3284.904 

488, 

.359 

3038, 

.602 

602.050 

3987.742 

581, 

.156 

3811, 

.123 

638.380 

4298.508 

628, 

.588 

4214, 

►  746 

638.993 

4303.754 
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.840 

4293, 

►  889 

604.556 

4009.174 

611, 

.644 

4069, 

►  808 

551.875 

3563.929 

563. 

.784 

3664, 

►  470 

488.422 
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503, 

.495 

3162, 

►  859 

430.014 

2565.149 

444, 

.704 

2682, 

►  313 

378.825 

2162.272 

392, 

.244 

2265, 

►  878 

334.335 

1822.364 

346, 

.389 

1912, 

►  664 

292.033 

1508.545 

303, 

►  778 

1593, 

►  458 

253.818 

1237.030 

264, 

,697 

1314, 

►  324 

221.745 

1020.835 

231, 

.493 

1086, 

.063 

192.665 

830.846 

202, 

.093 

889, 

.326 

162.528 

648.094 

172, 

.390 

707, 

.893 

136.630 

499.075 

146, 

.123 

550, 

.939 

112.799 

372.226 

122, 

.143 

419, 

,922 

88.473 

256.190 

98, 

.883 

303, 

.215 

72.484 

188.034 

80, 

.855 

221, 

.774 

65.718 

161.139 

70, 

►  075 

178, 

.460 

63.653 

152.931 

65, 

►  344 

159, 

►  655 

63.085 

150.673 

63, 

.648 

152, 

.914 

62.949 

150.134 

63, 

.116 

150, 

.796 

62.915 

150.000 

62, 

.962 

150, 

►  185 

62.915 

150.000 

62, 

.915 

150. 

.000 

10 
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REACH    STORAGE  -  OUTFLOW   RELATIONSHIP 
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Figure  3, — Reach  storage-outflow  relationship. 
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Table  3. — Input  for  example  2 


3600  3600  50   50 

1   27 

35 

001 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

00 

40 

50 

60 

80 

100 

120 

150 

200 

250 

2 

00 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

00 

800 

900 

1000 

1100 

1200 

1300 

1400 

1500 

4 

ion   1 

0 

.1633 

.5179 

1.6369 

4.1364 

8.7595 

16.2298 

29,43'S4 

45.4092 

1 

11 

72.3008 

103.49 

138.41 

217.91 

308.26 

408.21 

572.  12 

875.32 

1209.9 

2 

11 

1565.15 

1939.71 

2325.76 

2724.55 

3133.9 

3548.1 

3970.2 

4397,9 

4831.0 

3 

11 

5709.5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

11 

0 

250 

650 

1188 

1775 

2775 

3775 

4275 

4400 

1 

22 

4275 

3838 

3338 

2775 

2338 

1963 

1650 

1338 

1088 

2 

22 

900 

713 

525 

400 

27 

0 

0 

0 

0 

3 

22 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-r05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

00 

40 

50 

60 

80 

100 

120 

150 

200 

250 

2 

on 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

on 

800 

900 

1000 

1100 

1200 

1300 

1400 

1500 

4 

ion    1 

0 

.1633 

.5179 

1.6369 

4.1364 

8.7595 

16.2298 

29.4354 

45.4092 

1 

11 

72.3008 

103.49 

138.41 

217,91 

308.26 

408.21 

572.  12 

875.32 

120^.9 

2 

11 

1565.15 

1939.71 

2325.76 

2724.55 

3133.9 

3548.1 

3970.2 

43'57,9 

48''1  .0 

-3 

1  1 

S709.5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

11 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

-.05 

7200 

0 

' 

002 

Table  '^. ---Input  for  example  shown  in  figure  5,  at  top 


3600  3600  50   50 

I   27 

35 

001 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

00 

40 

50 

60 

80 

100 

120 

150 

200 

250 

2 

00 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

on 

800 

900 

1000 

1100 

1200 

1300 

1400 

1500 

4 

100    1 

0 

.1633 

.5179 

1.6369 

4.1364 

8.7595 

16.2298 

29.4354 

45.4092 

1 

11 

72.3008 

103.49 

138.41 

217.91 

308.26 

408.21 

572.12 

875.32 

1209.9 

2 

11 

1565.15 

1939.71 

2325.76 

2724.55 

3133.9 

3548.1 

3970.2 

4397.9 

4831.0 

3 

11 

5709.5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

11 

0 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

150 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

22 

0- 

0 

0 

0 

0 

0 

0 

0 

0 

2 

22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3 

22 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

00 

40 

50 

60 

80 

100 

120 

150 

200 

250 

2 

00 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

on 

800 

900 

1000 

HOC 

1200 

1300 

1400 

1500 

4 

100    1 

0 

•  1633 

.  .5179 

1.6369 

4.1364 

8.7595 

16.2298 

29.4354 

45.4092 

1 

n 

72.3008 

103. *9 

1)8.41 

217.91 

308.26 

408.21 

572.12 

875.32 

1209.9 

2 

1  ] 

1565.15 

1939,71 

2325.76 

2724.55 

3133.9 

3548.1 

3970.2 

4397.9 

48^1. 0 

•s 

11 

3709.5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

11 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3 

22 

7200 

0 

002 
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problem.  Inflow  is  a  uniform  and  constant  lateral  rate  of  0.005  c.f.s./ft.  The  rate  could  be 
assumed  to  be  falling  over  a  watershed  and  could  be  varied  so  as  to  evaluate  the  influence  of 
temporal  and  spatial  variations  of  rainfall.  Program  input  is  presented  in  table  5. 

Table  5. — Input  for  example  shown  in  figure  5,   at  bottom 


3600  3600  50   50 

1   27 

35 

001 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

1  00 

40 

50 

60 

8C 

100 

120 

150 

200 

250 

2 

1  00 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

1  00 

800 

900 

1000 

1100 

1200 

1300 

1400 

1500 

4 

100. 

0 

.1633 

.5179 

1.6369 

4.1364 

8.7595 

16,2298 

29.4354 

45,4092 

1 

1  11 

72.3008 

103.49 

138.41 

217.91 

308,26 

408.21 

572.12 

875. ?2 

1209,9 

2 

1  n 

1565.15 

19(39.71 

2325.76 

2724.55 

3133.9 

3548,1 

3970.2 

4397.9 

4831.0 

3 

1  11 

5709.5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

1  11 

.0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1  22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

1  22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3 

1  22 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

0 

0 

1 

2 

4 

7 

11 

16 

23 

30 

1 

1  00 

40 

50 

60 

80 

100 

120 

150 

200 

250 

2 

1  00 

300 

350 

400 

450 

500 

550 

600 

650 

700 

3 

1  0-) 

800 

900 

1000 

1100 

1200 

1300 

1400 

1500 

4 

100 

0 

•  1633 

.5179 

1.6369 

4.1364 

8.7  595 

16,2298 

29.-4354 

45.4092 

1 

]  11 

72.3008 

103.49 

138.41 

217.91 

308.26 

408.21 

572,12 

875,32 

1209.9 

2 

1  ]  1 

1565.15 

1939,71 

2325.76 

2724.55 

3133.9 

3548.1 

3970.2 

4397,9 

4831 .0 

3 

1  n 

5  70  9'.  5 

6601 

7503 

8408 

9326 

10245 

11177 

12103 

4 

]  11 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

•  005 

.005 

.005 

.005 

,005 

.005 

,005 

.005 

.005 

•  005 

.005 

.005 

.005 

.005 

.005 

0 

7200 

0 

002 

CONSISTENCY  OF  THE  STORAGE  FLOOD  ROUTING  METHOD 

A  vexing  characteristic  of  storage  flood  routing  methods  is  the  dependence  of  calculated 
discharges  on  the  selection  of  reach  lengths  (Ax)  and  time  increments  (At).  An  indication  of 
this  difficulty  with  respect  to  storage  flood  routing  with  coefficients  Muskingum  was  discussed 
by  Gilcrest  3,  thus:  "However,  it  has  been  found  that  actual  flood  hydrographs  along  several 
hundred  miles  of  river  can,  through  judicious  selection  of  routing  reaches  and  time  interval  used 
in  routing,  be  duplicated  quite  accurately. ..." 

Several  numerical  trials  were  made  to  evaluate  the  influence  of  Ax  and  At  on  the  storage 
flood  routing  without  coefficients  method.  Tables  6  and  7  represent  the  routing  of  the  same  inflow 
hydrograph  used  in  example  1.  Three  reach  lengths  are  combined  with  two  time  increments.  The 
indication  is  that  the  effect  of  Ax  over  the  range  shown  here  for  a  constant  At  is  small.  However, 
an  influence  of  At  is  clearly  evident. 

Table  8  presents  a  set  of  routings  for  a  triangular- shaped  hydrograph  inflow  into  a  prismatic 
channel.  The  conclusion  is  that  the  peak  reduction  is  completely  dependent  on  At. 

This  conclusion  was  first  pointed  out  by  Thomas  (9).  Especially  evident  from  table  8  is  a 
possible  trend  toward  no  crest  reduction.  This  is  in  agreement  with  the  work  of  Henderson  (4) 
on   kinematic   waves.    Thus,    the   method   of   storage  flood  routing  is  converging  on  the  correct 

result.    In   view   of   the  approach  to  a  condition  of  no  peak  reduction  as  At ►  0,  an  effort  is 

required    to   select  At   so   as   to   approximate   the  true  peak  reduction  condition.  Presented  in 


14 


Table  6. --Comparison  of  hydrographs  for  a  time  increment  of  1  hour 
with  three  reach  lengths 


Time 

Hydrograph  discharge 

( hoiirs ) 

Ax  =  3,600  feet 

4x  =  7,200   feet 

Ax   =  14,400  feet 

0 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

22 

24 

26 

28 

30 

32 

34 

36 


C.f.s. 

150 

183 

560 

1,866 

3,502 

4,117 

3,665 

2,804 

2,038 

1,443 

1,001 

667 

419 

263 

187 

158 

151 

150 

150 


C.f.s. 

150 

182 

560 

1,866 

3,503 

4,117 

3,665 

2,804 

2,038 

1,443 

1,001 

667 

419 

263 

187 

159 

151 

150 

150 


C.f.s. 

150 

179 

561 

1,867 

3,507 

4,118 

3,663 

2,802 

2,038 

1,443 

1,001 

667 

420 

261 

187 

159 

152 

150 

150 


Table  7. --Comparison  of  hydrographs  for  a  time  increment  of  one-half  hour 

with  three  reach  lengths 


Time 

Hydrograph  discharge 

(hovirs) 

Ax   =  3,600  feet 

Ax  =  7,200  feet 

Ax  =  14,400   feet 

0 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

22 

24 

26 

28 

30 

32 

34 

36 


C.f.s. 

150 

158 

443 

1,768 

3,622 

4,262 

3,788 

2,816 

2,023 

1,435 

994 

662 

410 

243 

169 

150 

150 

150 

150 


C.f.s. 

150 

157 

445 

1,769 

3,623 

4,262 

3,767 

2,816 

2,024 

1,435 

994 

662 

411 

243 

169 

152 

150 

150 

150 


C.f.s. 

150 

150 

452 

1,771 

3,628 

4,^60 

3,765 

2,814 

2,024 

1,435 

995 

662 

412 

241 

171 

153 

150 

150 

150 
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Table  8. --Comparison  of  hydrographs  for  several  time  intervals  with  a 
constant  reach  length  (prismatic  channel) 


Time 

Inflow 

Outflow  72,000  feet  downstream  for  4t  of-- 

(hoijrs) 

1/10  hour 

1/4  hour 

1/2  hour 

1  hour 

C.f.s, 


C.f.s. 


C.f.s. 


C.f.s, 


C.f.s, 


0 

150 

150 

150 

150 

150 

1/2 

600 

150, 

150 

150 

-- 

1 

1,200 

150 

150 

150 

159 

1  1/2 

1,800 

150 

150 

151 

-- 

2 

2,400 

150 

150 

159 

260 

2  1/2 

3,000 

150 

150 

204 

-- 

3 

3,600 

150 

177 

391 

733 

3  1/2 

4,200 

286 

531 

827 

-- 

4 

4,800 

1,276 

1,330 

1,476 

1 

,721 

4  1/2 

4,390 

2,260 

2,206 

2,207 

-- 

5 

3,975 

3,045 

2,989 

2,898 

2 

,801 

5  1/2 

3,575 

3,741 

3,636 

3,451 

-- 

6 

3,180 

4,277 

4,060 

3,801 

3 

,460 

6  1/2 

2,790 

4,402 

4,193 

3.933 

-- 

7 

2,400 

4,194 

4,075 

3,875 

3 

,551 

tables  9  and  10  are  the  comparisons  of  the  selection  of  At  in  the  storage  flood  routing  procedure 
so  as  to  correspond  to  the  peak  reduction  found  with  a  higher  order  routing  method  (1).  A  At  in 
the  range  of  one-fourth  to  one-half  hour  appears  to  give  a  satisfactory  correspondence.  Figures 
6  and  7  show  the  tabular  comparison  between  routed  hydrographs  for  two  different  hydrograph 
shapes.  In  figure  7,  the  1/4- hour  hydrograph  could  not  be  shown  since  it  fell  too  close  to  the 
hydrograph  predicted  by  the  higher  order  routing  method. 


CONCLUSIONS 

Formulation  of  the  unsteady,  non-uniform  open  channel  flow  problem  in  pure  kinematics 
leads  to  a  simple  problem  for  solution.  The  computer  program  listed  here  appears  to  be  func- 
tioning quite  well  over  a  wide  spectrum  of  problems. 

An  unavoidable  feature  of  storage  flood  routing  methods  is  still  of  concern;  namely,  the 
selection  of  the  reach  length  (Ax)  and  time  increment  (At).  Tentative  results  of  the  numerical 
testing  in  this  report  indicate  that  the  only  concern  in  selecting  the  reach  length  is  to  fit  the 
physical  dimension  of  the  channel  system.  However,  the  time  increment  selection  has  major 
influence  on  the  shape  of  the  routed  hydrograph.  For  two  extreme  inflow  hydrograph  shapes, 
several  time  increment  comparisons  were  made  with  a  higher  order  routing  method,  i.e.,  a  more 
complete  equation  of  motion  was  utilized.  The  best  correspondence  occurred  when  the  time 
increment  was  At    =  tp/16 


where 


tp    =    time  to  peak  of  the  inflow  graph. 


Future  development  for  this  routing  method  will  require  extensive  testing  with  data  from 
actual  floodflow  situations.  It  may  happen  that  the  strongly  nonprismatic  character  of  actual 
channels  may  reduce  the  influence  of  time  increments. 


16 


Table    9. --Comparison  of  storage   flood  routing  read   in   1/2-hour  time    increment 
without    coefficients  with  a   higher   order   formulation 


Time 

Inflow 

Outflow' 

Outflow^ 

( hours) 

4 1  =  1/4  hour 

At  =   1/2  hour 

C.f.s. 


C.f.s. 


0 

150 

150 

1/2 

600 

150 

1 

1,200 

150 

1  1/2 

1,800 

150 

2 

2,400 

150 

2  1/2 

3,000 

150 

3 

3,600 

159 

3  1/2 

4,200 

390 

4 

4.800 

1,380 

4  1/2 

4,390 

2,455 

5 

3,980 

3,126 

5  1/2 

3,580 

3,703 

6 

3,180 

4,002 

6  1/2 

2,790 

4,062 

7 

2,400 

3,923 

C.f.s. 

150 

150 

150 

150 

150 

150 

177 

531 

1,330 

2,206 

2,989 

3,636 

4,060 

4,193 

4,075 


C.f.s. 

150 

150 

150 

151 

159 

204 

391 

827 

1,476 

2,207 

2,898 

3,451 

3,801 

3.933 

3,875 


'  Calculated  with  the  equation  of  motion,  neglecting  only  the  inertia  forces. 
^  Calculated  by  storage  flood  routing  without  coefficients. 


Table  10, --Comparison  of  storage  flood  routing  read  in  hourly  time  increments 
without  coefficients  with  a  higher  order  formulation 


Time 
( hours ; 


Inflow 


Outflow' 


Outflow 


At  sz   1/4  hour 


At   =   1/2  hour 


At   =  1   hour 


C,f.s. 


C.f.s, 


C.f.s, 


C.f.s, 


C,f,s, 


0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 


150 

150 

150 

150 

150 

400 

154 

150 

150 

150 

650 

157 

150 

157 

182 

1,075 

182 

150 

213 

288 

1,775 

349 

380 

445 

559 

2,775 

921 

910 

949 

1,082 

3,775 

1,810 

1,730 

1,769 

1,866 

4,275 

2,837 

2,790 

2,767 

2,753 

4,400 

3,737 

3,700 

3,623 

3,503 

4,275 

4,176 

4,190 

4,119 

3,970 

3,875 

4.274 

4,320 

4.262 

4.117 

3,375 

4,101 

4,150 

4,119 

3,988 

2,775 

3,742 

3,800 

3,767 

3,665 

Calculated  with  the  equation  of  motion,  neglecting  only  the  inertia  forces. 
Calculated  by  storage  flood  routing  without  coefficients. 
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Symbol 
Input: 

DELT 

DELTI 

DELA 

DELAl  - 

TOLR 

Nl 

N2 

AREAl(I) 
DISCHl(I) 
AREA2(I) 
DISCH2(I) 


SYMBOLS  USED  IN  THE  PROGRAM 

Description 

Time  increment  used  in  routing  computations  (constant)  in  seconds 

Same  as  DELT,  saved  for  initialization 

Area  increment  in  rating  table  (see  Solve  subroutine) 

Same  as  DELA,  saved  for  initialization 

Tolerance  (see  Solve  subroutine) 

Number  of  inflow  hydrograph  entries 

Number  of  tabulated  rating  function  entries 

Entries  in  the  upper  and  lower  section  rating  function  tables 
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Symbol 

(I     =     1,  N2) 

QUI) 

QLKI) 

QL2(I) 

(I     =    1,  Nl) 

Q2(I) 

(I     =    1) 

DELX 

Output: 

LLl 

LL2 

Aid) 

QUI) 

A2(I) 

Q2(I) 

(I     =    1,  M) 

Main  Program: 

READ 

TBLLP 

XINT 

(CD, 

E,  F,  G) 

Description 


Entries  in  the  inflow  hydrograph 


Entries  in  the  upper  and  lower  lateral  inflow  hydrographs 


Initial  flow  rate  at  the  lower  section 

Distance  increment  (can  be  changed  for  each  reach) 

Upper  section  number 
Lower  section  number 

Upper  section  flow  area  and  flow  rate 
Lower  section  flow  area  and  flow  rate 


Nil 

M 

I 


See  Read  subroutine 

See  Table  Look-up  subroutine 

Function  for  interpolation 
Arithmetic  statement  arguments 

C  -  Lower  bound  for  wanted  table 
D  -  Upper  bound  for  wanted  table 
E  -  Lower  bound  for  argument  table 
F  -  Upper  bound  for  argument  table 
G  -  Argument 

Value  of  Nl  saved  for  initialization 
Total  number  of  outflow  entries 
Index  of  do  loop 
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Symbol 
Jl 

J 

JC 

KKl  (1  or  2) 

ALPHA 

BETA 

XI 

SOLVE 

N3 

Read  Subroutines: 

A(I1) 

(II    =    1,  Kl) 

II 

Kl 


Description 

Lower  index  value  of  a  do  loop  used  for  filling  out  lateral  inflow  hydrograph 
entries 

Index  code  for  incrementing  I 

Index  code  for  decrementing  I 

Index  code  for  a  computed  go  to 

First  term  of  the  right-hand  side  of  equation  6 

Second  term  of  the  right-hand  side  of  equation  6 

Sum  of  ALPHA  and  BETA 

See  Solve  subroutine 

Number  of  flow  entries  occurring  after  inflow  ceases 


List  items  of  entries  being  read 

Index  of  implied  do  loop 
Number  of  items  in  list 


Table  Look- Up  Subroutine: 


X4 

A  (J) 

(J    =   1,  N2) 

J 

K 

Solve  Subroutine; 
A2(I) 
Q2(I) 
X2 
X2 
AU 
FAU 
AL 
FAL 
DELA 
ANEW 
X3 
TOLR 


Look-up  argument 
Look-up  table 

Upper  index  of  table  bracket 
Lower  index  of  table  bracket 

Current  area  value  for  iteration 

Current  flow  value  for  iteration 

Left-hand  side  of  equation  6 

Difference  between  X2  and  XI 

Upper  bound  of  area  in  iteration  process 

X2  if  value  of  area  too  large 

Lower  bound  of  area  in  iteration  process 

X2  if  value  of  area  too  small 

Area  increment  or  decrement  used  to  establish  an  upper  or  lower  bound 

New  value  of  area  in  iteration  process 

Difference  between  new  area  value  and  previous  area  value 

Iteration  cutoff 
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APPENDIX  1 
MAIN  LINE  PROGRAM 


c 

STORAGE  FLOOD  ROUTING  WITHOUT  COEFFICIENTS 

oni 

c 

»»♦»«»»*»»»»#*#»#»•#»»♦#*#*»»*»#*»#»*»♦♦*»*#*##»»##»»#»*♦#»»)<■»♦*#» 

30  2 

DIMENSION  AREA  1(50)  .DISCHK  50  )  »Q  1  (  1  00  )  .  Al(  lOU  )  .  AREA.^  I  bU  ) 

no3 

DIMENSION  DISCH2(50) .02 t 100 ) . A2 (  100  )  .OL  111 UO ) .QL2 ( 10  0 ) 

004 

COMMON  AREAl.DISCHl.Ol.Nl.I . J » K .M. AREA2 »DI SCH2 . A2 .02 . X 1 .DEL  1 .DEL! 1 

nn5 

COMMON  N2.X4.X2»DELA.DELAl.ANEi^.TOLR.X3 

006 

COMMON  A1.DELX.QL1.QL2 

007 

c 

«♦♦#♦#»♦♦##»«»»»♦»»♦»»•»#*#«»»»♦♦»♦«#»•###»♦♦##»#»**#♦»»»**»**»#»* 

008 

XINT(C»D.E.F.G)=C+( (G-E)/(F-E) ) » ( D-C ) 

UU9 

c 

C=LOWER  BOUND  FOR  WANTED  TABLE 

ulu 

c 

D=UPPER  BOUND  FOR  WANTED  TABLE 

on 

c 

E=LOWER  BOUND  FOR  ARGUMENT  TABLE 

012 

c 

F=UPPER  BOUND  FOR  ARGUMENT  TABLE 

013 

c 

G=ARGUMENT 

014 

c 

*♦♦♦»»#♦»»*##»♦»♦♦»♦»#♦♦»»»»»**»*♦♦#»♦»»»♦•##»♦*»♦»♦#♦*»»»»»•»**»• 

015 

c 

READ  PAR  CARD 

016 

READ  l.DELT.DELT I.DEL A. DELAl.TOLR.Nl .N2 

017 

1 

FORMAT(2F5.0.3FA.1.2I3) 

018 

c 

»»»•♦♦«»»•#«»#««««»♦♦#»»»#»*#♦*»»#•**♦**♦»*♦*»*»»»»»#»*»♦*♦♦♦♦♦#*» 

019 

c 

READ  FIRST  SECTION  RATING  AND  INFLOW 

02n 

CALL  READ  IAREA1.N2) 

021 

CALL  READ  <DISCH1,N2) 

022 

CALL  READ  (01. Nl) 

023 

CALL  READ(OLl.Nl) 

-^24 

c 

»♦«♦«♦♦«»♦•»»»#»»♦#»#»»»♦»**##*»♦#*##**»•#«#»♦*♦»»«*•*»»*»#♦**♦♦»* 

025 

c 

CALCULATE  FIRST  SECTION  AREAS 

026 

19 

DO  10  1=1. Nl 

027 

CALL  TBLLP(DISCHl.Ql) 

028 

10 

A1<I)«XINT(AREA1((C).AREA1(J).DISCH1(K).DISCH1(J)»Q1(I)) 

029 

c 

»»»»#»#»»»#*♦»#»»♦*♦*♦#»»»»*#»*****♦»»**»»***»»*«»***»»*»**♦«»♦»*• 

330 

c 

READ  SECOND  SECTION  RATING  AND  LATERAL  INFLOW 

331 

LL1  =  0 

032 

LL2  =  0 

033 

N11=N1 

034 

59 

M  =  N1 

035 

99 

1  =  1 

336 

N1=N11 

037 

J1=N1+1 

038 

IF(M-Nl)  952.952.951 

039 

951 

DO  950  J=J1»M 

040 

950 

QL2( J)=0.0 

041 

952 

CONTINUE 

042 

CALL  READ  (AREA2.N2) 

043 

CALL  READ(DISCH2.N2) 

044 

CALL  READ(QL2.N1) 

345 

c 

#»»*»»»»»#»«•»»♦•»♦♦»»»#♦*«♦»»»♦•****♦*♦»**»•*♦*♦*»»**•»###»»»**»* 

046 

c 

READ  SECOND  SECTION  INITIAL  VALUE  AND  DELX 

047 

READ  2.DELX»Q2( I ) 

048 

2 

FORMAT  (F5.0.F8.0) 

049 

CALL  TBLLP(DISCH2tQ2) 

050 

A2( I »«XINT(AREA2(K) .AREA2(J) .DISCH2(K) .DISCH2( J) .Q2( I  )  ) 

051 

J=l 

052 

N1«M 

053 

c 

•***«*»«•**«*•*»••••**•«****««*******««*«****************»»««««*«« 

054 

c 

ROUTING  DURING  INFLOW 

055 

100 

ALPHA«(  AK  I  )+A2(  1  )  )/2. 

056 

I=«J+1 

057 

BETA«(DELT/DELX}*Q1(I )♦( -Al ( I )  +  ( DELT ) » ( QLl ( I)+0L2(  I)  )  )/2. 

058 
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APPENDIX  1— Continued 
MAIN  LINE  PROGRAM— Continued 


X1=ALPHA+BETA  059 

IF  (XI) 101.101 tlOZ  ■  060 

101  Q2( I )=Q2(1 )  061 

A2( I  )=A2(  1)  062 

DELT=DELT+DELT1  063 

I=I-J  064 

J=J+1  ,  065 

IF(J-N1 )1C0, 105*105  066 

105  1=J  067 

GO  TO  111  .     /        '  06  8 

102' CALL  SOLVE  ]  069 

IF  (Q2( I)-Q2(l)-. 1)698. 698. 699  070 

698  Q2(  I )=Q2(  1)  071 
A2(I)»A2(1)  072 

699  CONTINUE  073 
IF(SENSE  SWITCH  2)  700.701  074 

700  PUNCH  7C2.A2« I ) .Q2< I ) .1  075 
702  FORMAT  (2F13,3.I3)  076 

701  CONTINUE  077 
DELT«0ELT1  •  .  .  078 
IF  ( I-N1)110.111.111  079 

110  J=I  080 
GO  TO  100  081 

Q  ••««•«••••••••••••••••**«***««*«***»**»«***««*««*«*»•*«*««««*«««««  0  8  2 

L  ROUTING  AFTER  INFLOW  ^^    -■  383 

111  N3»0  0  84 

112  J=I  085 
N3«N3-fl  086 

501  ALPHA»(Al(Nn+A2(  I  )  )/2.  087 
I«J+1  088 
BETA=(DELT/DELX)*01(N1)-K-A1(N1  )+DELT*(QLl(Nl  )+QL2(Nl  )  )  )/2.  089 
Xl»ALPHA*bETA  090 
IF(X1)502»502.505                                                   .r  091 

502  Q2( I )«02(1)  092 
A2(I)»A2(1)  093 
GO  TO  600  094 

505  CALL  SOLVE  095 

IF  (02(  n-Q2n)-. 1)502. 502. 601  ^'  096 

601  GO  TO  112  097 

600  M=M+N3  098 

IF  (M-100)602.602.801  099 

801  M=100  100 

602  CONTINUE  101 
I2-N1+1  102 
00  900  I«I2»M  103 
QLK  I)«0.0  104 
0L2( 1 )»0.0  105 
A1(I)«A1(N1)  106 

900  on  I  )«01(N1)  1''7 

LL1«LL1*1  i  108 

LL2»LL1*1  119 

C      »««*»««««««««««««««••»«»»««#*«•«*««*•«••»«*«»«««««««**«•»«•«*«•«««  110 

C      PUNCH  OUT  AND  INTERCHANGE  111 

PUNCH  4  112 

4  FORMAT(///)  113 
PUNCH  5.LL1.LL2  H'' 

5  F0RMAT(16H  IN  SECTION  NO  =I3.14X.17H  OuT  SECTION  NO  =13)  }15 
PUNCH  6  116 
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APPENDIX  l~Continued 
MAIN  LINE  PROGRAM— Continued 


6    F0RMAH4X.8H     IN    AREA»6X.9H     IN    DISCH,10X.9H    OUT    ARFA.5X.10H    OUT    DIS  117 

iCH)  lie 

3  FORMAT  (2F13,3»5X.2F13.3)  11<? 

DO  115  I=1»M  120 

PUNCH  3tAl ( I ) *Q1( I ) tA2( I ) tOZi I )  121 

QLHI)»QL2(n  122 

A1(I)=A2(I»  123 

115  OK  I  )«Q2(  I  )  124 

C      INTERCHANGE  RATING  TABLES  126 

DO  116  I«1.N2  127 

AREAKI  J«AREA2«  I  )  128 

116  DISCHK  I  )=0ISCH2(  n  129 
GO  TO  99  130 

800  END  131 


Read  and  Table  Look-up  Subroutine 


READ  SUBROUTINE 

SUBROUTINE  READ(A*K1) 

DIMENSION  A( 100) 

DIMENSION  AREA H 50). DISCHK 50) .01 ( 100 ) . Al( 100 ) .AREA2(50) 

DIMENSION  DISCH2(S0)*O2( 100) •A2 ( 100 ) >QL1( 100 ) *QL2 ( 100 ) 

COMMON  AREAl»DISCHl»01»Nl*I*J*K*M»AREA2»DISCH2tA2«Q2«XltDELI *DELI 1 

COMMON  N2tX4«X2»0ELA»DELAl»ANEW»T0LR.X3 

COMMON  Al»OELX*OLltOL2 

READ  2»(A(I1)*I1«1»K1) 

FORMAT  (9F8.0) 

RETURN 

END 


C      TABLE  LOOK-UP  SUBROUTINE 
C      A-ARGUMENT  TABLE 
C      B«ARGUMENT 

SUBROUTINE  TBLLPCA.B) 

DIMENSION  A( 100) .B(50) 

DIMENSION  AREA  1(50) .DISCHK  50) .01 (100 ) . Al  (  inc ) » AREA2 ( 50 ) 

DIMENSION  DISCH2<  50 ) .02 ( 1 00 ) .A2 ( IOC ) .QL K 1 00 ) .012 ( 100 ) 

COMMON  AREA1.0ISCH1.Q1»N1.I.J.K»M.ARFA2,DI SCH2 » A2 , 02 . X 1 .DFLT.DrLTl 

COMMON  N2.X4.X2.DELA.0ELA1,ANEW»T0LR.X3 

COMMON  A1.DELX.QL1.QL2 

X4  =  B(  I  ) 

IF  (X4)16.15.16 

15  J  =  2 

GO  TO  20 

16  DO  30  J=1.N2 

IF  (A(J)-X4)30.20.20 
20  K=J-1 

RETURN 
30  CONTINUE 

END 


001 

002 

003 

004 

005 

0C6 

007 

008 

009 

010 

on 

012 

001 

002 

003 

004 

005 

006 

007 

008 

009 

010 

Oil 

012 

013 

014 

015 

016 

017 

018 

019 

02  0 
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APPENDIX  l—Continued 
Solve  Subroutine 


c 

SOLVE  ROUTINE 

001 

SUBROUTINE  SOLVE 

002 

DIMENSION  AREA1(50).DISCH1(50) .0 1( 100  )  » Al ( 100 ) .AREA2 ( 50 ) 

003 

DIMENSION  D1SCH2(50) .02 ( 100 ) . A2 ( 100 ) .QL 1 ( 100) .0L2 ( 100 ) 

004 

COMMON  AREA1.DISCH1,Q1,N1»I  »  J.K.  .M,AREA2  .D  I  SCH2  ♦A2  .02  tX  1  »DELT  .DELT  1 

005 

COMMON  N2»X4.X2»DELA»DELA1,ANEW.T0LR.X3 

006 

COMMON  A1.DELX.QL1.QL2 

007 

XlNT(C»D.E.F,G)«C-»-(  (G-E)/(F-E)  )»(D-C) 

008 

AU  =0.                                        —:'.'■ 

009 

FAU  =0. 

010 

f 

AL  =  0. 

on 

FAL  =0. 

012 

DEL  A  «DELA1 

013 

IF(I-Nl)  800.800»801 

014 

801 

A2(I)«A1(N1) 

015 

GO  TO  11 

016 

800 

A2( I  )*A1(  I  ) 

017 

11 

CALL  TBLLPJAREA2.A2) 

018 

02 ( I )»X1NT(DISCH2(K).DISCH2( J)»AREA2(K),AREA2(J).A2( I ) ) 

019 

IF(SENSE  SWITCH  1)700.701 

020 

700 

PUNCH  702.A2(n.O2(  I). I 

021 

702 

FORMAT  (2F13.3.I3) 

022 

701 

CONTINUE 

023 

X2«(DELT/DELX)«Q2n  )  +  <A2(  I  )  )/2. 

024 

X2«X2-X1                                               ^ 

025 

IF  (X2)3.5.2 

026 

2 

DELA=DELA1 

027 

AU«A2( I )                                                    .  ' 

028 

FAU=X2 

029 

IF  (AL)4.12«4                                                       ':.    , 

;   030 

12 

A2(I »=A2(I)-DELA                                       '        .  . 

031 

13 

IF  <A2( I ))7.7.11 

032 

7 

DELA=DELA».5 

033 

A2( 1 )«A2( I )*OELA 

034 

GO  TO  13 

035 

3 

DELA=DELA1 

036 

AL=A2{I > 

037 

FAL=-X2 

038 

IF  (AU)4.14.4 

039 

14 

A2( I )xA2(I)+DELA 

040 

15 

IF  (A2m-AREA2«N2)  »  11»11  •16 

041 

16 

DELAxDELA».5 

042 

A2( I )«A2( I )-DELA 

043 

GO  TO  15 

044 

4 

ANEW»AU-(FAU/(FAU*FAL» )*(AU-AL) 

045 

X3«ANEW-A2<n 

046 

X3=ABSF(X3) 

047 

IF  (X3-TOLR)5.5.10 

048 

10 

A2(I)«ANEW 

049 

GO  TO  11 

050 

5 

A2(I)«ANEW 

051 

CALL  TBLLP(AREA2.A2) 

052 

02  (  I  )«XINT(DISCH2(IC)»DISCH2(  J).AREA2(K)  .AREA2(  J)  »A2{  I  )  ) 

053 

RETURN 

054 

END 

055 
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APPENDIX  2 


Program  Input  Format' 


Parameter  Card  1: 


1-5                         6-10              11-14              15-18 
XXXXX           XXXXX           XXX.X           XXX.X 
DELT              DELTI             DELA             DELAl 

19-22 
XXX.X 
TOLR 

23-25             26-28 
XXX               XXX 
Nl                   N2 

Hydrograph  and  Rating  Table  Entries: 

1-8                      9-16 
XXXXXXXX.      XXXXXXXX. 

17-24 
XXXXXXXX. 

25-32 
XXXXXXXX, 

33-40                    41-48 
XXXXXXXX.       XXXXXXXX. 

49-56                        57-64 
XXXXXXXX.       XXXXXXXX. 

65-72 
XXXXXXXX. 

73-80^ 
XXXXXXXX. 

Parameter  Card  2; 

1-5                  6-13 
XXXXX.         XXXXXXXX. 

Inflow  and  Outflow  Values: 


Program  Output  Format""" 


1-13 

XXXXXXXXXX.XXX 

Al(I) 


14-26 
XXXXXXXXXX.XXX 
Ql(I) 


32-44 
XXXXXXXXXX.XXX 
A2(I) 


45-57 
XXXXXXXXXX.XXX 
Q2(I) 


Loading  Sequence 

1.  Parameter  card  1 

2.  Flow  area  portion  of  inflow  section  rating  table 

3.  Flow  portion  of  inflow  section  rating  table 

4.  Inflow  hydrographs  at  inflow  section 

5.  Lateral  inflow  at  inflow  section 

6.  Flow  area  portion  of  outflow  section  rating  table 

7.  Flow  portion  of  outflow  section  rating  table 

8.  Lateral  inflow  at  outflow  section 

9.  Parameter  card  2 


^Decimal  is  shown  at  its  implied  position.  If  another  position  is  required,  it  should  be  punched, 
^Columns  73-80  are  reserved  for  identification,  such  as  card  nuniber,  card  type,  section  number,  watershed 
location,  etc. 
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